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Serotonin (5-HT) is implicated in the development of stress-related mood disorders 
in humans. Physical activity reduces the risk of developing stress-related mood 
disorders, such as depression and anxiety. In rats, 6 weeks of wheel running protects 
against stress-induced behaviors thought to resemble symptoms of human anxiety and 
depression. The mechanisms by which exercise confers protection against stress-induced 
behaviors, however, remain unknown. One way by which exercise could generate stress 
resistance is by producing plastic changes in gene expression in the dorsal raphe 
nucleus (DRN). The DRN has a high concentration of 5-HT neurons and is implicated 
in stress-related mood disorders. The goal of the current experiment was to identify 
changes in the expression of genes that could be novel targets of exercise-induced 
stress resistance in the DRN. Adult, male F344 rats were allowed voluntary access to 
running wheels for 6 weeks; exposed to inescapable stress or no stress; and sacrificed 
immediately and 2 h after stressor termination. Laser capture micro dissection selectively 
sampled the DRN. mRNA expression was measured using the whole genome Affymetrix 
microarray. Comprehensive data analyses of gene expression included differential gene 
expression, log fold change (LFC) contrast analyses with False Discovery Rate correction, 
KEGG and Wiki Web Gestalt pathway enrichment analyses, and Weighted Gene 
Correlational Network Analysis (WGCNA). Our results suggest that physically active rats 
exposed to stress modulate expression of twice the number of genes, and display a more 
rapid and strongly coordinated response, than sedentary rats. Bioinformatics analyses 
revealed several potential targets of stress resistance including genes that are related 
to immune processes, tryptophan metabolism, and circadian/diurnal rhythms. 

Keywords: Affymetrix gene microarray. Weighted Gene Correlational Network Analysis, bioinformatics, laser 
capture microdissection, stress resistance, dorsal raphe nucleus 



INTRODUCTION 

Depression and anxiety frequently coexist and are the most 
common mood disorders affecting society. The World Health 
Organization estimates that 121 million people currently suffer 
from depression. Individuals suffering from depression have sig- 
nificant impairment in quality of life (Rapaport et al., 2005), are 
at increased risk for developing coronary heart disease (Wulsin 
and Singal, 2003) and type 2 diabetes (Knol et al., 2006), and have 
higher mortality due to suicide. By 2030, depression is expected 
to be a leading cause in the global burden of disease (Mathers and 
Loncar, 2006). 

Stressful life events often precede the onset of depression 
(Kendler et al, 1999; van Praag, 2005) and anxiety. Despite the 
high occurrence and significant disability associated with stress- 
related mood disorders, the pathophysiology of these conditions 
is not fully understood. Important to note is that not every 



individual who experiences a stressful life event develops a seri- 
ous mood disorder, and these individuals may possess resistance 
to the negative affective consequences of stress. Pinpointing the 
factors by which stress resistance occurs could provide a bet- 
ter understanding of the neurobiological mechanisms underlying 
stress-related mood disruptions. 

To investigate the neural circuitry underlying stress-related 
mood disorders, researchers use animal models (Krishnan and 
Nestler, 2008). Rats exposed to an acute inescapable stressor, such 
as tail shock, later exhibit behaviors argued to resemble symptoms 
of human anxiety and depression (Maier and Watkins, 1998), and 
these behaviors are responsive to pharmaceutical treatment with 
anxiolytics (Drugan et al., 1984) and antidepressants (Sherman 
et al., 1982). Inescapable stressor exposure also produces various 
physiological perturbations. Long-term increases in basal levels 
of plasma corticosterone and decreases in corticosteroid-binding 
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globulin occur in rats following tail shock (Fleshner et al., 1995). 
Additionally, acute stress increases interleukin-ip (IL-ip), leading 
to immune modulation (Moraska et al., 2002), and centrally, con- 
tributes to behavioral consequences of stress (Maier and Watkins, 
1995). Circadian-regulated processes are also susceptible to acute 
stress. Thompson et al. (2013) observed a decrease in amplitude 
and disruption in diurnal pattern of core body temperature and 
heart rate in rats exposed to tail shock. Moreover, inescapable 
stress produces alterations in brain serotonergic circuits. The 
serotonergic system has long been implicated in underlying the 
behavioral consequences of inescapable stress exposure in rats 
(Maier and Watkins, 2005) and has been heavily implicated in 
human affective disorders (Sharp and Cowen, 2011). 

Numerous components of the serotonergic system such as 
serotonin (5-HT) receptors, the 5-HT transporter, and extra- 
cellular 5-HT levels are sensitive to stress. Serotonergic nerve 
terminals and receptors also occupy regions of the brain involved 
in neuroendocrine and behavioral responses to stress (Chaouloff, 
1993). One region of particular interest is the dorsal raphe 
nucleus (DRN), a small midbrain structure containing a high 
concentration of stress-responsive 5-HT cell bodies (Grahn et al., 
1999). Hyper activation and sensitization of DRN 5-HT neurons 
is thought to underlie the depression- and anxiety-like behav- 
iors induced by inescapable stress exposure (Maier et al., 1995; 
Christianson et al, 2008). 

The DRN receives afferent, and provides efferent, projections 
to brain regions involved in fear, anxiety, and depression. These 
regions include the prefrontal cortex, striatum, bed nucleus of 
the stria terminalis (BNST), amygdala, and locus coeruleus (LC). 
Efferent DRN projections render these regions susceptible to 
stress-induced 5-HT activity in the DRN. Furthermore, these 
regions are themselves sensitive to stress (Cullinan et al., 1995), 
provide afferent input to the DRN, and may modulate DRN 5-HT 
activity. Nerve terminals containing corticotropin-releasing factor 
(CRF), a neuropeptide produced in response to elevated Cortisol 
levels, for example, are present in the DRN (Swanson et al., 1983). 
Given that the BNST projects to the DRN and contains many CRF 
neurons (Day et al., 1999), the BNST is believed to be a primary 
source of CRF to the DRN. Interestingly, CRF injected into the 
DRN increases 5-HT activity in a subpopulation of cells (Lowry 
et al., 2000), and injection of CRF into the caudal DRN produces 
behaviors resembling those produced by inescapable stress expo- 
sure (Hammack et al, 2002). Thus other brain regions influence 
DRN 5-HT levels, and interactions between those regions and the 
DRN likely contribute to the DRN's role in stress-related mood 
disorders. 

Also important to consider is that within the DRN, inter- 
actions between diverse cell populations may influence stress- 
induced 5-HT activity. The DRN is not just a homogenous 
structure of 5-HT neurons. Other populations of neurons con- 
taining the neurotransmitters y-aminobutyric acid (Belin et al., 
1979; Day et al, 2004), dopamine (Lindvall and Bjorklund, 1974; 
Stratford and Wirtshafter, 1990), and glutamate (Commons et al., 
2005) also exist. Cells containing neuropeptides such as sub- 
stance P (Hokfelt et al, 1978) and neuropeptide Y (de Quidt 
and Emson, 1986) are also present. These various neuropeptides 
and neurotransmitters/receptors are capable of modulating 5-HT 



(Ferre et al., 1994; Song et al, 1996; Tao and Auerbach, 2000; 
Valentino et al, 2003). Therefore, stress-induced alterations in 5- 
HT activity within the DRN and at DRN projection sites may be 
influenced indirectly through non-serotonergic neuronal modu- 
lation of serotonergic neurons. Non-serotonergic neurons in the 
DRN are also sensitive to 5-HT, and can have inhibitory and exci- 
tatory responses to 5-HT release (Marinelli et al, 2004). Dynamic 
interactions between serotonergic and non-serotonergic neurons 
originating at DRN afferent sites and within the DRN likely con- 
tribute to the effect of stress on net DRN 5-HT release within the 
DRN and at DRN projections sites. 

Non-neuronal cell types, such as astrocytes and microglia, 
may also influence DRN neural activity. Microglia are the res- 
ident "immune cells" of the brain and are sensitive to stress- 
induced elevation of glucocorticoids (Nair and Bonneau, 2006; 
Sugama et al, 2007). Activated microglia release interleukin- 1 
(IL-1) (Giulian et al., 1986), tumor necrosis factor-a (TNF-a) 
(Sawada et al, 1989), and interleukin-6 (IL-6) (Righi et al, 1989). 
Inescapable stress increases IL-ip in the brain (Nguyen et al., 
1998). Stress-induced activation of microglia may occur in the 
DRN and effect 5-HT neurons. Consistent with this idea, admin- 
istration of interferon-y (IFN-y) and TNF-a reduced the survival 
of 5-HT neurons in organotypic DRN sections (Hochstrasser 
etal.,2011). 

Overall, the DRN is an important region of investigation in 
studying the neurobiological mechanisms of stress-related mood 
disorders. Elucidation of variables influencing the serotonergic 
response to stress within the DRN may provide a better under- 
standing of the development of these disorders. Furthermore, 
identification of interventions that prevent or manipulate the 
serotonergic response to stress and/or influence the various fac- 
tors capable of modulating 5-HT activity within the DRN, may 
lead to the identification of novel therapeutic targets. 

In humans, physical activity is one factor known to influ- 
ence an individual's response to stress. Exercise reduces the risk 
of developing stress-related depression and anxiety (Fox, 1999). 
Similarly, in rats, 6 weeks of voluntary wheel running protects 
against the behavioral consequences of inescapable stress expo- 
sure (Greenwood et al., 2003). It is believed that wheel running 
prevents these behaviors by attenuating stress-induced activation 
of 5-HT neurons within the DRN. Wheel running may do this 
by producing plasticity at (1) DRN afferent sites (2) DRN effer- 
ent sites or (3) within the DRN itself (Greenwood and Fleshner, 
2011). Given that hyperactivity of 5-HT neurons in the DRN 
is necessary for the development of stress-induced behaviors 
in rats and our lab has previously shown that wheel running 
attenuates stress-induced c-fos expression in DRN 5-HT neu- 
rons (Greenwood et al., 2003), we will focus on exercise-induced 
plastic changes that may occur within the DRN itself. 

In particular, the 5-HTi^ inhibitory autoreceptor has been 
implicated in the mechanism by which wheel running could 
constrain stress-induced 5-HT activity and protect against the 
behavioral consequences of inescapable stress. 5-HTia receptors 
inhibit the activity of 5-HT neurons (Sprouse and Aghajanian, 
1987) and reduce 5-HT release (Casanovas et al., 1997). Six 
weeks of wheel running increases 5-HTi^ mRNA expression in 
the DRN (Greenwood et al., 2003, 2005) and thus, may increase 
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5-HTi^ receptor-mediated inhibition of DRN 5-HT neurons 
during inescapable stress. 

The protective effect of wheel running could also occur 
indirectly, through a non-serotonergic route. One possibility is 
through neuropeptides. Wheel running increases brain-derived 
neurotropic factor (BDNF) (Neeper et al., 1995), a neuropep- 
tide important for maintaining neuronal health and function, 
and galanin (Tong et al., 2001) in the hippocampus, and also 
upregulates gene expression of galanin in the LC (Holmes et al., 
2006; Sciolino et al., 2012). Wheel running may also increase 
levels of BDNF and galanin in the DRN. Both factors are coex- 
pressed in 5-HT neurons in the DRN (Merlio et al., 1992; Xu 
and Hokfelt, 1997) and are capable of modulating 5-HT activ- 
ity. Infusion of BDNF into the DRN modifies the neuronal firing 
of 5-HT by decreasing the regularity of the firing pattern (Celada 
et al., 1996). Additionally, an in vitro study revealed that galanin 
hyperpolarizes 5-HT neurons within the DRN (Xu et al., 1998). 
Exercise-induced increases in BDNF and galanin may protect 
against stress-induced activation of 5-HT neurons through mod- 
ulating and, in the case of BDNF, inhibiting 5-HT neuronal 
activity. 

Another method by which exercise may confer protection is 
through an immune-related mechanism. Evidence suggests a role 
of cytokines in human mood disorders (Maes, 2008) and stress- 
induced behaviors in rats. Injection of an IL-1 receptor antagonist 
into the brain blocks stress-induced depression- and anxiety- 
like behaviors in rats (Maier and Watkins, 1995), suggesting that 
activity at brain IL-1 receptors is important for the production 
of these behaviors. Speaker et al. (2011) observed that 6 weeks 
of wheel running attenuates stress-induced increases in plasma 
IL-ip, one of two cytokines that bind IL-1 receptors. It is pos- 
sible that 6 weeks of wheel running also reduces stress-induced 
increases in brain IL-ip, and through reducing ligand availability, 
protects against stress-induced alterations in brain IL-1 recep- 
tor activity. Given that the DRN contains many IL-1 receptors 
(Cunningham and De Souza, 1993), it may be particularly sensi- 
tive to stress-induced and/or exercise-induced alterations in IL-1 
receptor activity. 

Though the precise mechanisms are not fully understood, 
the protective effect of exercise likely involves preventing stress- 
induced alterations in the serotonergic system, either by directly 
constraining activity of 5-HT neurons within the DRN or indi- 
rectly, through altering other neurotransmitter systems or neu- 
ropeptides within the DRN that are capable of modulating 5-HT 
neurons. Furthermore, DRN 5-HT neurons may be influenced 
by exercise-induced plastic changes that reduce afferent input 
to the DRN, activate afferent inhibition of the DRN during 
stress (Greenwood and Fleshner, 2011), or produce alterations in 
postsynaptic 5-HT receptor function (Greenwood et al., 2012). 
Elucidating the mechanism by which exercise produces stress- 
resistance and protects against the behavioral consequences of 
stress may lead to the identification of novel therapeutic targets 
and development of more targeted drugs for the treatment of 
human stress-related mood disorders. 

One approach to reveal novel targets is by employing the 
use of microarray technology. Microarray technology permits 
the investigation of the expression of tens of thousands of genes 



simultaneously, at the level of mRNA transcription. Predesigned 
chips that contain sequences, known as probes, derived from 
every gene within a specified genome can be probed with mRNAs 
obtained from experimental samples in order to gain informa- 
tion about gene expression under the given conditions (Cox et al., 
2012). When used in conjunction with laser capture microdissec- 
tion, microarrays can reveal expression patterns of genes within 
specific cells. Using microarray and laser capture microdissection, 
therefore, it is possible to assess the effect of exercise and/or stress 
on gene expression in cell populations specific to the DRN. 

The purpose of this experiment was to investigate the effect 
of exercise and/or stress on gene expression within the DRN. 
We hypothesized that wheel running produces changes in mRNA 
transcription within the DRN, and physically active rats exposed 
to stress have different gene expression profiles compared to 
sedentary rats exposed to stress. The differences in gene expres- 
sion patterns within the DRN between physically active and 
sedentary rats exposed to stress may underlie the molecular 
mechanisms by which exercise protects against behaviors pro- 
duced by inescapable stress exposure. Whole genome Affymetrix 
microarray analysis was used to assess gene expression. Our goal 
was to use an exploratory approach to ( 1 ) systematically organize 
the transcriptome (17,170 genes) obtained from the microar- 
ray analysis into a more manageable and focused gene set and 
(2) extrapolate physiological implications from this focused gene 
set by identifying novel targets of exercise-induced stress resis- 
tance within the DRN. To ensure a comprehensive assessment of 
the data, the organizational process involved two approaches, (1) 
identification of genes based on changes in differential expres- 
sion in response to exercise and/or stress (2) identification of 
genes based on changes in coexpression in response to exer- 
cise and/or stress. For the differential expression analysis, two 
measures of significance were utilized. In a more conservative 
approach, genes were identified by log fold changes in gene 
expression. In a second, less stringent approach, genes statistically 
significantly differentially expressed by p < 0.05 were identified. 
These p-values were corrected for multiple comparisons using the 
False Discovery Rate adjustment method. The coexpression anal- 
ysis narrowed the transcriptome from 17,170 genes to 1 1 modules 
of highly coexpressed networks of genes. These networks of genes 
were then correlated to the response to exercise and/or stress. 
Both the differential and coexpression analysis returned sets of 
genes that were further sorted by their relationship to functional 
categories derived from bioinformatics databases. Novel targets 
of exercise-induced stress resistance were identified within these 
functional categories. 

MATERIALS AND METHODS 
ANIMALS 

The University of Colorado Boulder Animal Care and Use 
Committee approved all protocols for this study. A total of 
48 adult, male Fisher 344 rats weighing 170-180 grams at 
time of arrival (Harlan Laboratories) were used in this experi- 
ment. Upon arrival, animals were individually housed in Nalgene 
Plexiglas cages (45 x 25.2 x 14.7 cm). The housing environment 
was maintained on a 12:12 h lightdark cycle, controlled for 
humidity, and held at a constant temperature of 22° C. Rats were 
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allowed ad libitum access to food and water and were weighed 
weekly to ensure each animal remained healthy. Following arrival, 
animals were acclimated to the housing conditions for 1 week 
before experimental manipulation. 

WHEEL RUNNING 

Animals were randomly assigned to remain sedentary (Sed, 
« = 23) or were housed with a running wheel (Run, n = 25), 
and allowed voluntary access to the wheel for 6 weeks. During 
the 1-week acclimation period, wheels were rendered immobile 
with metal stakes. Daily wheel revolutions each animal ran were 
logged with Vital View software (Mini Mitter). The product of 
the total number of daily revolutions and the wheel circumfer- 
ence (1.081 m) was calculated to obtain daily running distance. 
Daily running distance was summed in order to get an average 
weekly running distance. 

INESCAPABLE STRESS 

Animals were randomly assigned to remain in their home cages 
(HC) or receive inescapable tail shock (Stress). The stress pro- 
cedure occurred between 0700 and 1200. Animals subjected to 
stress were restrained in acrylic cylinders (23.4 x 7 cm diameter). 
The tail projected from the back of the restraint device. An elec- 
trode was positioned 3 cm from the base of the tail and served as 
the vehicle by which shock was delivered. The shock procedure 
consisted of 100, 5-s tail shocks administered on a random 60-s 
inter- trial interval. Rats received 1.0 mA tail shocks for 50min, 
at which time the intensity of shock was increased to 1.5 mA 
tail shocks for the remainder of the session. The entire stress 
procedure lasted 1 h and 48 min. Rats were sacrificed by rapid 
decapitation immediately following termination of tail shock 
(StressO) or 2 h post termination of tail shock (Stress2). The sacri- 
ficing of rats that remained in their home cage was time matched 
with those animals subjected to tail shock. 

TISSUE COLLECTION AND CRY0SECTI0NING 

RNAse free conditions were maintained throughout tissue pro- 
cessing. After rats were sacrificed, brains were extracted and 
flash frozen at — 20°C, in 2-Methylbutane (Fisher Scientific), for 
4 min. Brains were stored at — 80°C prior to sectioning. Brains 
were prepared with M-l embedding matrix before sectioning 
at — 21°C with a cryostat (Leica CM1850). Tissue was sectioned to 
a thickness of 20 um through the rostral to mid-caudal (approx- 
imately —7.3 to —8.2 mm relative to Bregma) portions of the 
DRN. This specific region of the DRN was targeted because 
it is involved in modulating stress- and anxiety-like behaviors 
(Hale et al., 2012) and prior evidence suggests that alterations 
in gene expression occur in this region following 6 weeks of 
wheel running (Greenwood et al., 2003, 2005). Sections were 
freeze-mounted to PEN membrane frame slides (MDS Analytical 
Technologies) and stored at — 80°C until further use. 

LASER CAPTURE MICRODISSECTION AND RNA ISOLATION 

Laser capture microdissection was performed to procure a pre- 
cise sample of the DRN from each rat. Slides containing sections 
of DRN were allowed to thaw for 20 s prior to being fixed in 
75% ethanol, subjected to a Histogene Stain (for visualization 
purposes), and dehydrated in graded ethanol concentrations, in 



accordance with the Arcturus Histogene LCM Frozen Section 
Staining Kit protocols (Applied Biosystems). Following staining 
procedures, slides were loaded into the laser capture microdissec- 
tion system (Arcturus XT, Life Technologies). The regions of DRN 
targeted for capture were the dorsal and ventral portions of the 
rostral to mid-caudal DRN. Samples were captured so that each 
sample contained the entire portion of the dorsal and ventral por- 
tion of the DRN at the given rostral-caudal level. DRN samples 
were obtained by using an infrared laser to adhere the tissue to a 
cap coated with a thermoplastic film (Capsure Macro LCM Caps, 
Applied Biosystems). An ultraviolet laser was used to separate the 
DRN from the rest of the brain section. An average of 23 DRN 
samples, ranging in size from 300,000 to 800,000 um 2 (depending 
on rostral to caudal level), were successfully dissected and pooled 
for each rat to ensure maximal total RNA yield. Following laser 
capture microdissection, caps were incubated in RNA extraction 
buffer (Applied Biosystems) for 30 min and frozen at — 80°C until 
future use. RNA was isolated using the Arcturus Picopure RNA 
Isolation Kit (Applied Biosystems) in accordance with kit proto- 
cols. Samples were stored in Elution Buffer (Applied Biosystems) 
at — 80° C until microarray analysis. 

MICR0ARRAY ANALYSIS 

Samples were sent to the Genomics and Microarray Core 
Facility at the University of Colorado Denver for whole genome 
analysis using microarray. RNA integrity was evaluated with 
the Agilent Bioanalyzer 2100 and RNA 6000 Nano/Pico Kit 
(Agilent Technologies). Concentrations of extracted RNA were 
assessed with the Nanodrop spectrophotometer (Nanodrop 
Technologies). One sample was removed from further processing 
due to poor integrity of RNA (n = 47). A total of 100-150 ug RNA 
per each sample was converted to double stranded cDNA and 
then transcribed into cRNA using the Ambion WT Expression 
Kit, in accordance with kit protocols. Following generation of 
cRNA, second cycle, first strand cDNA synthesis was carried out 
in order to transform the cRNA into single-strand cDNA. The 
cDNA was fragmented and the Genechip WT Terminal Labeling 
Kit (Affymetrix) was used to label the single-stranded DNA with 
biotin. Samples were hybridized to an Affymetrix Genechip Rat 
Gene 1.1 ST Array Platform. Hybridization, washing, staining, 
and scanning were executed using the GeneTitan instrument 
(Affymetrix). 

MICROARRAY DATA PRE-PROCESSING 

The Bioconductor toolset within the R statistical software pro- 
gram was used to format the raw microarray data. This pre- 
processing was completed using the 'expresso' option in the 'affy' 
package of the Bioconductor toolset and included background 
adjustment, log fold transformation, and normalization. To con- 
trol for inter-array variability, the dataset was normalized using 
the Robust Multi Array Average method. Gene chip and RNA 
quality were assessed by examining total mRNA expression for 
each animal. 

MICROARRAY CONTRAST GENERATION 

Following pre-processing and normalization, a standardized 
expression value was obtained for each gene for each rat. 
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The expression values for each gene were averaged for each 
experimental group. The LIMMA package was used to gen- 
erate nine contrasts between experimental groups. These con- 
trasts included [(Sed S t re ssO v. Sed H c) v. (Run Stl -essO v. Run H c)], 
[(Sed Stre ss2 v. Sed H c) v. (Runs tre ss2 v. Run H c)]> Run H c v. Sed H c, 

SedstressO V. Sed H C, Sed St ress2 v. Sed H c, Run St ressO v. Run H c, 
Run St ress2 V. Run H C, Run S tressO V. SedsressO, and Run St ress2 v. 

Sedsress2- For each contrast, the difference in the expression level 
of each individual gene was calculated by subtracting the expres- 
sion level of the 2 nd group in each contrast from the expression 
level of the 1 st group in each contrast. For example, the con- 
trast Runnc v. SedHc indicates that the expression level of gene 
X in the Sednc group was subtracted from the expression level of 
gene X in the Runnc group, or (RunHc — Sednc)- The first two 
contrasts represent the interaction between exercise and stress at 
each time point. For each contrast, p-values and test statistics were 
calculated for each gene according to the absolute value of differ- 
ence in gene expression observed between the groups. The False 
Discovery Rate multiple-test adjustment method was applied in 
the calculation of these p-values in order to control for the chance 
of yielding false positive (significant) results. The log fold change 
(LFC) in gene expression was also calculated for each gene in each 
contrast. Out of 27,000 possible genes, 17,170 gene transcripts 
were reliably detected. These genes were considered the transcrip- 
tome, or the genes expressed in cells of the DRN as a result of the 
experimental conditions. 

DIFFERENTIAL GENE EXPRESSION IDENTIFICATION AND 
BIOINFORMATICS 

In an initial approach, genes differentially expressed by a LFC 
> ±1.1 were identified for each contrast. A second approach was 
performed utilizing the same contrasts as previously described. 
However, less stringent requirements for statistical significance 
were utilized (p < 0.05) to identify differentially expressed genes 
between groups. Genes that were significantly differentially 
expressed at a p < 0.05 were organized into nine sets, one set 
for each contrast, and imported into the bioinformatics system, 
Web Gestaldt. Specifically, KEGG (Kanehisa and Goto, 2000) and 
Wikipathways (Wild) (Kelder et al., 2012) pathway enrichment 
analysis were applied to each gene set in order to identify the top 
functionally enriched pathway categories related to the genes sig- 
nificantly differentially expressed in each contrast. Both KEGG 
and Wiki databases were used in an effort to generate a more com- 
prehensive analysis. The KEGG system is recognized as one of the 
major pathway databases (Bauer-Mehren et al, 2009), whereby 
data is derived from published work. KEGG pathway content 
includes categories in metabolism, genetic information process- 
ing, organismal systems, environmental information processing, 
cellular processes, and drug development. The Wiki system is 
curated by the scientific community and serves as a complemen- 
tary and enhancing source to the KEGG database (Bauer-Mehren 
et al., 2009). Finally, ANOVAs were performed on select genes of 
interest that were identified through the pathway analysis. 

WEIGHTED GENE CORRELATIONAL NETWORK ANALYSIS 

Given that genes often operate in a coordinated manner to 
accomplish a physiological function, a more sophisticated 



approach utilizing Weighted Gene Correlational Network 
Analysis (WGCNA) was also performed. That is, in the absence 
of absolute differences in gene expression, the coexpression of 
genes may differ across experimental conditions. The WGCNA 
package within the R software program was used to perform this 
analysis. Following standard preprocessing and normalization 
of the data, a gene expression profile was available for each rat. 
Based on this expression profile, rats were clustered hierarchically 
within a dendogram based on Euclidian distance, or similarity 
between expression profiles. The dendogram was visualized to see 
how the physical traits (experimental conditions) related to the 
various clusters. Next, modules of highly coexpressed genes were 
identified and related to physical traits. Importantly, the genes 
within each module are more highly correlated with each other 
than to the rest of the transcriptome. Physical traits were catego- 
rized by experimental group (Sed.HC, Sed.StressO, Sed.Stress2, 
Run.HC, Run.StressO, Run.Stress2) and differences between 
groups (RunvsSed.HC, StressOvsHC.Sed, Stress2vsHC.Sed, 
StressOvsHC.Run, Stress2vsHC.Run, RunvsSed.StressO, 
RunvsSed.Stress2). A correlation value and p-value associ- 
ated with the strength of correlation was calculated for each 
module. These values are considered a representation of the 
correlation and correlational strength of the module to each 
physical trait. Modules with a correlational strength of p < 0.001 
were targeted for further investigation by ANOVA. Following 
ANOVA analysis, modules that had statistically significant main 
effects of exercise, stress and/or an exercise by stress interaction 
were subjected to KEGG and Wiki analyses. 

RESULTS 

BODY WEIGHT AND RUNNING DISTANCE 

A repeated measures ANOVA was used to analyze body weights. 
Repeated measures ANOVA analysis revealed statistically signif- 
icant main effects of time [F(6, 252) = 518.415; p < 0.0001] and 
exercise [F(i : 42) = 7.759; p = 0.0080] and a reliable time by exer- 
cise interaction [F^ 252) = 2.634; p = 0.0171] on body weight. 
Running distance steadily increased over the course of the exper- 
iment from 5959.305 ± 382.081 m during week 1 to 19355.603 ± 
2983.808 m during week 6 [% U 5) = 18.870; p < 0.0001]. 

ASSESSMENT OF RNA AND MICROARRAY GENECHIP QUALITY 

To verify microarray chip quality and mRNA integrity, boxplots 
were constructed that represented total mRNA expression for 
each rat. Visual inspection of boxplots revealed four outliers. One 
sample was previously identified with spectrophotometry anal- 
ysis. The additional three outliers were dropped from further 
analysis. Final group totals were Sednc {n = 6), SedstressO (n = 7), 

Sed S tress2 (« = 7), Run H C (« = 8), Run St ress0 (« = 8), Run Str ess2 

(n = 8), for a total of (« = 44) rats. 

DIFFERENTIAL GENE EXPRESSION ANALYSIS RESULTS 

The effect of wheel running and/or exposure to stress on log fold 

changes in gene expression of ±1.1 in the DRN 

In a conservative initial approach, genes differentially expressed 
by a LFC > ± 1 . 1 in response to exercise and/or stress were identi- 
fied. Overall, relatively few genes had LFCs in expression > ±1.1. 
The effect of stress on LFCs in gene expression > ±1.1 is not 
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different between sedentary and physically active rats immedi- 
ately following or 2h after stress exposure. When considering 
the effect of exercise, only one gene was statistically significantly 
altered. This gene was transthyretin, which was downregulated 
in physically active rats. Following stress exposure, transythretin 
remained downregulated in physically active rats immediately 
after, but not 2 h post stress. Compared to home cage controls, 
stress produced alterations in gene expression regardless of phys- 
ical activity status at both time points. Figure 1 shows the effect 
of stress on total number of genes altered (A) and in which direc- 
tion (B) by a LFC > ±1.1 in sedentary rats and physically active 
rats immediately following and 2 h post stress relative to home 
cage non-stressed controls. Exposure to stress produced alter- 
ations in gene expression immediately following and 2 h after in 
both sedentary and physically active rats. A greater number of 
genes changed in physically active rats (41 and 27) than seden- 
tary rats (18 and 26) for both time points (Figure 1 A). All 18 
genes that were differentially expressed immediately following 
stress in sedentary rats were also differentially expressed immedi- 
ately following stress in physically active rats. Twenty-two of the 
twenty-six genes that were differentially expressed 2 h post stress 
in sedentary rats were also differentially expressed 2 h post stress 
in physically active rats. Two genes were differentially expressed 



exclusively in the physically active rats exposed to stress. These 
genes were CD180 and fos-like antigen 1, which were different 
at both time points. There were also two genes that were differ- 
entially expressed exclusively in sedentary rats exposed to stress 
compared to home cage non-stressed controls. These genes were 
only differentially expressed 2 h post stress and were fibronectin 1 
and solute carrier organic anion transporter family, member lcl. 

Figure IB shows the number of genes that were upregulated 
and downregulated by a LFC > 1.1 in sedentary and physically 
active rats in response to stress. For all groups, a greater num- 
ber of genes were differentially upregulated than downregulated 
in stressed rats compared to home cage non-stressed controls. 

The effect of wheel running and/or exposure to stress on changes 
in gene expression ofp< 0.05 in the DRN 

A second less stringent approach utilized p-values that were 
corrected for multiple comparisons to identify genes that were 
differentially expressed bjp < 0.05 in response to exercise and/or 
stress. The effect of stress on differential gene expression is differ- 
ent depending on physical activity status. Differential expression 
of 1028 genes was observed immediately following stress and 
637 genes were differentially expressed 2 h post stress in seden- 
tary rats compared to physically active rats. (These results are 
detailed in Figure 3 and Table 3). Compared to sedentary rats, 
physically active rats had differential expression in 2350 genes 
(1290 upregulated, 1060 downregulated). Immediately follow- 
ing stress, differential expression of 634 genes (413 upregulated, 
221 downregulated) was observed in physically active rats com- 
pared to sedentary rats. Two hours post stress, 997 genes (610 
upregulated, 387 downregulated) were differentially expressed in 
physically active rats compared to sedentary rats. Compared to 
home cage non-stressed controls, stress produced changes in gene 
expression in both sedentary and physically active rats. Figure 2 
shows the effect of stress on total number of genes altered (A) 
and in which direction (B) by p < 0.05 in sedentary rats and 
physically active rats immediately following and 2 h post stress 
relative to home cage non-stressed controls. Exposure to stress 
produced alterations in gene expression immediately following 
and 2 h after stress in both sedentary and physically active rats. 
For both time points, physically active rats exposed to stress had a 
greater number of genes differentially expressed relative to home 
cage non-stressed controls than sedentary rats exposed to stress 
(Figure 2A). For both sedentary and physically active rats, stress 
altered the expression of a greater number of genes 2 h after com- 
pared to immediately following stress. This difference was greater 
in sedentary rats. 

Figure 2B indicates the number of genes that were upregulated 
and downregulated by p < 0.05 in sedentary and physically active 
rats in response to stress. For both sedentary and physically active 
rats, a greater number of genes were differentially upregulated 2 h 
after compared to immediately following stress. In contrast, the 
number of genes differentially downregulated was greater imme- 
diately following stress compared to 2 h after for both sedentary 
and physically active rats. For both time points, physically active 
rats had a greater number of genes upregulated than sedentary 
rats. Similarly, physically active rats had a greater number of genes 
downregulated than sedentary rats at both time points. 
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FIGURE 1 | Number of genes differentially expressed by a log fold 
change > ±1.1 in rats exposed to stress compared to home cage 
non-stressed controls. (A) Total number of genes differentially expressed 
and (B) Number of genes differentially upregulated or downregulated. 
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FIGURE 2 | Number of genes differentially expressed by p < 0.05 in 
rats exposed to stress compared to home cage non-stressed controls. 

(A) Total number of genes differentially expressed. (B) Number of genes 
differentially upregulated or downregulated. 



KEGG functionally enriched pathway categories related to genes 
differentially expressed in the DRN following 6 weeks of wheel 
running and/or exposure to stress 

Genes that were differentially expressed at p < 0.05 across the 
various contrasts were imported into Web Gestaldt bioinformat- 
ics system and subjected to analysis with KEGG functional terms. 
Table 1 shows the top ten functionally enriched pathway cate- 
gories, for each contrast, related to genes differentially expressed 
following 6 weeks of wheel running and/or exposure to stress. 
Gene count refers to the number of genes from the data sets that 
contribute to each functional category. The p-value represents 
the statistical significance of each functionally enriched category 
identified. 

KEGG analysis revealed that genes that were differentially 
expressed between sedentary and physically active rats in response 
to stress were related to functional categories including metabolic 
pathways, mitogen-activated protein kinase (MAPK) signaling, 
neuroactive ligand receptor interaction, transforming growth 
factor-P (TGF-f5) signaling, epidermal growth factor family of 
receptor tyrosine kinases (ErbB) signaling, and vascular endothe- 
lial growth factor ( VEGF) signaling immediately following and or 
2 h post stress. 

Six weeks of wheel running modulated the expression of genes 
involved in physiological processes including metabolic activity, 



Table 1 | KEGG functionally enriched pathway categories generated 
from genes significantly differentially expressed at p < 0.05 in the 
DRN in response to exercise and/or stress. 



Contrast 

FUNCTIONALLY ENRICHED PATHWAY 



Gene count P-value 



(Sedstresso vs. Sed HC ) vs. (Run stress0 vs. Run HC ) 



Olfactory transduction:04740 


70 


1.65e-13 


Ribosome:03010 


18 


9.81 e-1 2 


Metabolic pathways:01100 


69 


1.48e-10 


MAPK signaling pathway:04010 


23 


6.84e-07 


Neuroactive ligand receptor interaction:04080 


20 


1.92e-05 


Pathways in cancer:05200 


23 


2.34e-05 


Endocytosis:04144 


17 


5.79e-05 


TGF beta signaling pathway:04350 


10 


6.16e-05 


ErbB signaling pathways:04012 


10 


0.0001 


Arachidonic acid metabolism:00590 


9 


0.0010 


(Sed str ess2 vs. Sed HC ) vs. (Run stress2 vs. Run 


He) 




Olfactory transduction:04740 


66 


1.78e-21 


Allograft rejection:05330 


8 


741e-06 


Pathways in cancer:05200 


18 


1.16e-05 


Autoimmune thyroid disease:05320 


8 


2.27e-05 


Graft-vs.-host disease:05332 


7 


5.39e-05 


C21-Steroid hormone metabolism:00140 


4 


4.44e-05 


Androgen and estrogen metabolism:00150 


6 


3.77e-05 


Type I diabetes mellitus:04940 


7 


0.0001 


Non-small cell lung cancer:05223 


6 


0.0003 


VEGF signaling pathway:04370 


7 


0.0003 


RunHc vs. Sed HC 


Metabolic pathways:01100 


141 


3.24e-17 


Olfactory transduction:04740 


115 


4.84e-12 


MAPK signaling pathway:04010 


48 


3.22e-12 


Pathways in cancer:05200 


53 


2.28e-11 


Ribosome:03010 


23 


6.46e-10 


FC epsilon Rl signaling pathway:04664 


21 


2.54e-09 


Cell cycle:04110 


27 


3.83e-09 


Long term depression:04730 


19 


5.84e-09 


Gap junction:04540 


21 


1.65e-08 


Vascular smooth muscle contraction :04270 


24 


4.09e-08 




MAPK signaling pathway:04010 


58 


1.10e-21 


Metabolic pathways:01100 


104 


2.33e-09 


VEGF signaling pathway:04370 


19 


5.06e-09 


Neurotropin signaling pathway:04722 


25 


9.26e-09 


Pathways in cancer:05200 


42 


2.58e-08 


Chronic myeloid leukemia:05220 


17 


4.77e-07 


Toll-like receptor signaling pathway:04620 


18 


7.50e-07 


GnRH signaling pathway:04912 


16 


1.93e-05 


Leukocyte transendothelial migration:04670 


18 


2.29e-05 


FC epsilon Rl signaling pathways:04664 


14 


2.75e-05 


Sed stress2 vs. Sed HC 


MAPK signaling pathway:04010 


56 


9.54e-19 


Pathways in cancer:05200 


53 


8.51e-13 


Adipocytokine signaling pathway:04920 


22 


1.95e-12 



(Continued) 
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Contrast 

FUNCTIONALLY ENRICHED PATHWAY 

Metabolic pathways:01100 
VEGF signaling pathway:04370 
Neuroactive ligand receptor interaction:04080 
Jak-STAT signaling pathway:04630 
Leukocyte transendothelial migration:04670 
Regulation of actin cytoskeleton:04810 
Toll-like receptor signaling pathway:04620 
Run s ,ressO vs. Run HC 
MAPK signaling pathway:04010 
Metabolic pathways:01100 
Pathways in cancer:05200 
Cytokine-cytokine receptor interaction:04060 
Adipocytokine signaling pathway:04920 
Focal adhesion:04510 
P53 signaling pathway:04115 
Neuroactive ligand receptor interaction:04080 
Calcium signaling pathway:04020 
Neutrotrophin signaling pathway:04722 
Run s , reS s2 vs. RunHc 
Pathways in cancer:05200 
MAPK signaling pathway:04010 
Metabolic pathways:01100 
Jak-STAT signaling pathway:04630 
Neuroactive ligand receptor interaction:04080 
Chronic myeloid leukemia:05220 
Cytokine-cytokine receptor interaction:04060 
Focal adhesion:04510 
Prostate cancer:05215 
Pancreatic cancer:05212 
RunstressO vs. Sed stress0 
Cytokine-cytokine receptor interaction:04060 
Pathways in cancer:05200 
Chemokine signaling pathway:04062 
MAPK signaling pathway:04010 
Toll-like receptor signaling pathway:04620 
Olfactory transduction:04740 
Apoptosis:04210 

Arachidonic acid metabolism:00590 
Jak-STAT signaling pathway:04630 
TGF beta signaling pathway:04350 
R"nstress2 vs. Sed stress2 

Neuroactive ligand receptor interaction:04080 

Prostate cancer:05215 

Pathways in cancer:05200 

Ribosome:03010 

Focal adhesion:04510 

Regulation of actin cytoskeleton:04810 

Melanoma:05218 

Olfactory transduction:04740 

Cell adhesion molecule:04514 

Wnt signaling pathway:04310 



Gene count P-value 



116 
19 

36 
25 
22 
31 
19 

63 

125 

50 

35 

17 

31 

17 

35 

28 

22 

74 

62 

127 

35 

44 

23 

37 

35 

23 

19 

12 

15 

9 

11 

6 

28 

6 

5 

7 

5 

23 
13 
24 
11 
16 
17 
9 

44 
12 
11 



3.16e-11 
1.68e-08 
1.80e-07 
2.22e-07 
3.05e-07 
4.89e-07 
4.78e-07 

2.05e-23 
9.13e-14 
5.84e-11 
5.59e-10 
6.23e-08 
8.65e-08 
3.64e-07 
7.63e-07 
1.25e-06 
3.15e-06 

3.25e-25 
3.30e-22 
8.15e-14 
1.15e-13 
4.00e-11 
6.70e-11 
6.25e-11 
8.43e-10 
1.01 e-09 
1.20e-08 

7.82e-05 

0.0003 

0.0018 

0.0038 

0.0038 

0.0032 

0.0047 

0.0058 

0.0086 

0.0103 

2.38e-07 

6.23e-07 

4.38e-06 

1.22e-05 

4.16e-05 

3.59e-05 

8.38e-05 

0.0004 

0.0006 

0.0011 



olfactory transduction, MAPK signaling, cell cycle, and long term 
depression. 

Compared to home cage non-stressed controls, both seden- 
tary and physically active rats exposed to stress had significant 
enrichment of functional categories related to MAPK signaling, 
metabolic pathways, adipocytokine signaling, and neuroactive 
ligand receptor interaction. Significant enrichment of functional 
categories including VEGF signaling and toll-like receptor sig- 
naling was exclusive to sedentary stressed rats compared to 
home cage non-stressed controls and occurred at both time 
points. Compared to home cage non-stressed controls, stress 
modulated the expression of genes involved in cytokine-cytokine 
receptor interaction exclusively in physically active rats at both 
time points. 

A direct comparison of sedentary and physically active rats 
exposed to stress, revealed enrichment differences in func- 
tional categories related to cytokine-cytokine receptor interac- 
tion, chemokine signaling, MAPK signaling, toll-like receptor 
signaling, apoptosis, janus kinase-signal transducer and activator 
of transcription (Jak-Stat) signaling, TGF-p 1 signaling, neuroac- 
tive ligand receptor interaction, cell adhesion molecules, and 
wingless-type mouse mammary tumor virus integration site 
(WNT) signaling either immediately following and/or 2h post 
stress. 

Wiki functionally enriched pathway categories related to genes 
differentially expressed in the DRN following 6 weeks of wheel 
running and/or exposure to stress 

Genes that were differentially expressed at p < 0.05 across the 
various contrasts were imported into Web Gestaldt bioinformat- 
ics system and subjected to analysis with Wiki functional terms. 
Table 2 shows the top ten functionally enriched pathway cate- 
gories, for each contrast, related to genes differentially expressed 
following 6 weeks of wheel running and/or exposure to stress. 
Gene count refers to the number of genes from the data sets that 
contribute to each functional category. The p-value represents the 
significance of each functionally enriched category identified. 

Wiki analysis revealed that genes that were differentially 
expressed between sedentary and physically active rats in response 
to stress were related to functional categories including metabolic 
pathways, MAPK signaling, adipogenesis, biosynthesis of aldos- 
terone and Cortisol, and diurnally regulated genes with circadian 
orthologs. In addition, various immune-associated categories 
were also identified including those related to the inflamma- 
tory and cytokine response as well as signaling pathways for 
interleukin-5 (IL-5), IL-6, B-cell receptor, TGF-f5 receptor, and 
TNF-a-nuclear factor kappa-light-chain-enhancer of activated 
B cells (NF-kB). 

Six weeks of wheel running modulated the expression of genes 
involved in physiological processes related to signaling path- 
ways for MAPK, epidermal growth factor receptor 1 (EGFR1), 
TNF-oi-NF-kB, Insulin, G-protein, IL-5, and B-cell receptor. 
Compared to home cage non-stressed controls, both sedentary 
and physically active rats exposed to stress had enrichment of 
functional categories related to signaling pathways for MAPK, 
insulin, TGF-P receptor, IL-6, EGFR1, delta notch, and toll-like 
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Table 2 | Wiki functionally enriched pathway categories generated 


Table 2 | Continued 






from genes significantly differentially expressed at p < 0.05 in the 








DRN in response to exercise and/or stress. 






Contrast 


Gene count P-value 


Contrast 


Gene count 


P-value 


FUNCTIONALLY ENRICHED PATHWAY 


















Toll-like receptor signaling 


16 


5.76e-07 


FUNCTIONALLY ENRICHED PATHWAY 






pathway:WP1309 






(Sed stress0 vs. Sed HC ) vs. (Run stress0 vs. Run 


He) 




IL-3 signaling pathway:WP319 


18 


5.39e-07 


Cytoplasmic ribosomal proteins:WP30 


17 


3.08e-09 


B cell receptor signaling pathway:WP285 


23 


1.21 e-06 


MAPK signaling pathway:WP358 


14 


2.17e-05 


Sedgtress2 vs. Sednc 






IL-5 signaling pathway:WP44 


9 


4.29e-05 


MAPK signaling pathway:WP358 


38 


6.26e-17 


Insulin signaling:WP439 


13 


7.33e-05 


Adipogenesis:WP1 55 


28 


5.99e-12 


B cell receptor signaling pathway:WP285 


13 


0.0001 


EGFR1 signaling pathway:WP5 


33 


1.42e-11 


Diurnally regulated genes with circadian 


6 


0.0005 


B cell receptor signaling pathway:WP285 


28 


r iQp_nQ 


orthologs:WP1306 






Insulin signaling:WP439 


27 


4 64p-09 


TGF beta receptor signaling 


11 


0.0005 


IL-3 signaling pathway:WP319 


21 




pathway:WP362 
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6 


0.0007 


pathway:WP1309 
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Dolta notch signaling pathway;WP199 


18 


5.10e-08 


IL-6 signaling pathway:WP135 


9 


0.0008 


TNF alpha NF-KB signaling pathway:WP457 


28 


6.69e-08 
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He) 
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0.0038 
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8.27e-13 
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0.0485 
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1 3 
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Kit receptor signaling pathway:WP147 


4 


0.0206 


orthologs:WP1306 






Inflammatory response pathway:WP40 


2 


0.0666 


Cardiovascular signalingiVVP590 


14 


3.15e-08 


Cytokines and inflammatory 


2 


0.0494 


Toll-like receptor signaling 


24 


7.80e-08 


response:WP271 






pathway :WP1 309 






Ovarian infertility genes:WP263 


2 


0.0535 


T cell rGCGptor signaling pathway:WP352 


23 


9.49e-08 


Metapathway biotransformation :WP1 286 


5 


0.0379 


Runstress2 V S- Runnc 






Run H c vs. Sed H c 






MAPK signaling pathway:WP358 


40 


4.75e-18 


MAPK signaling pathway:WP358 


29 


2.69e-09 


AdipogGnGsis:WP1 55 


33 


1.25e-15 


EGFR1 signaling pathway:WP5 


29 


3.80e-08 


EGFR1 signaling pathway:WP5 


39 


2.01 e-1 5 


TNF alpha NF-KB signaling pathway:WP457 


30 


3.17e-08 


Insulin signaling:WP439 


35 


2.93e-14 


Insulin signaling:WP439 


25 


4.32e-07 


B cell rGCGptor signaling pathway:WP285 


33 


5.96e-12 


Renin-angiotensin system :WP376 


13 


3.97e-07 


ILr6 signaling pathway;WP135 


26 


5.32e-12 


Myometrial relaxation and contraction 


24 


5.98e-07 


Delta notch signaling pathway:WP199 


21 


3.62e-10 


pathways:WP140 






TGF beta signaling pathways:WP505 


17 


6.00e-10 


Regulation of actin cytoskeleton:WP351 


24 


6.89e-07 


GPCRs, class A rhodopsin-like:WP473 


34 


5.30e-09 


G protein signaling pathways:WP73 


18 


1.14e-06 


TGF beta receptor signaling 


26 


6.09e-09 


IL-5 signaling pathway:WP44 


15 


3.31 e-06 


pathway:WP362 






B cell receptor signaling pathway:WP285 


24 


5.24e-06 


Rurictrpccn vs. Sedc+rp^n 






SedstressO vs. Sed H c 






Delta notch signaling pathway:WP199 


8 


2.30e-05 


MAPK signaling pathway:WP358 


36 


2.63e-16 


Kit receptor signaling pathway:WP147 


7 


6.36e-05 


Insulin signaling:WP439 


31 


1.25e-12 


IL5 signaling pathway:WP44 


7 


7.03e-05 


TGF beta receptor signaling 


25 


2.96e-09 


IL-3 signaling pathway:WP319 


7 


0.0007 


pathway:WP362 






IL-6 signaling pathway:WP135 


7 


0.0007 


GPCRs, class A rhodopsin-like:WP473 


32 


4.84e-09 


TGF beta signaling pathways:WP505 


5 


0.0011 


AdipogenesisiWPI 55 


23 


6.55e-09 


Toll-like receptor signaling 


6 


0.0012 


EGFR1 signaling pathway:WP5 


26 


767e-08 


pathway:WP1309 






IL-6 signaling pathway:WP135 


19 


1.53e-07 


Notch signaling pathway:WP517 


4 


0.0017 
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Table 2 | Continued 



Contrast 

FUNCTIONALLY ENRICHED PATHWAY 



Gene count 



P-value 



Endochondral ossification:WP1308 5 

Hedgehog signaling pathway:WP574 3 
R"nstress2 vs. Sed stress2 

Adipogenesis:WP155 11 

GPCRs, class A rhodopsin-like:WP473 1 5 

B cell receptor signaling pathway:WP285 12 

Hypothetical network for drug 5 
addiction:WP1281 

IL-6 signaling pathway:WP135 9 

Regulation of actin cytoskeleton:WP351 11 

Calcium regulation in the cardiac 10 
cell:WP326 

Cytoplasmic ribosomal proteins:WP30 9 

Androgen receptor signaling 9 
pathway:WP68 

Myometrial relaxation and contraction 10 
pathways:WP140 



0.0017 
0.0026 

0.0001 
0.0002 
0.0004 
0.0007 

0.0006 
0.0006 
0.0009 

0.0017 
0.0019 

0.0019 



Stress 0 



Stress 2 




FIGURE 3 | Effect of physical activity status on stress-induced 
differential gene expression in the DRN by p < 0.05. Number of genes 
differentially expressed immediately following (red) and 2 hours post (blue) 
stress, and genes differentially expressed at both time points (purple). 



receptor. Other genes affected by stress were related to func- 
tional categories including G protein coupled receptors (GPCRs) 
and adipogenesis. Compared to home cage non-stressed controls, 
stress modulated the expression of genes involved in interleukin-3 
(IL-3) signaling exclusively in sedentary rats at both time points. 
Significant enrichment of categories related to apoptosis, diur- 
nally regulated genes with circadian orthologs, cardiovascular 
signaling, and T-cell receptor signaling was exclusive to physically 
active stressed rats compared to home cage non-stressed controls. 

A direct comparison of sedentary and physically active rats 
exposed to stress, revealed significant differences in functional 
categories related to signaling pathways including delta notch, kit 
receptor, IL-3, IL-5, IL-6, TGF-f5 receptor, toll-like receptor, and 
B-cell receptor. Significant enrichment of functional categories 
related to adipogenesis and GPCRs was also observed. 



Table 3 | Functionally enriched KEGG and Wiki pathway categories 
related to genes differentially expressed in physically active and 
sedentary rats immediately following and 2 h post stress. 



Pathway analysis 


Gene count 


P-value 


FUNCTIONALLY ENRICHED PATHWAY 




KEGG ENRICHMENT ANALYSIS 


Olfactory transduction:04740 


19 


2.80en-08 


Pathways in cancer:05200 


8 


0.0001 


Metabolic pathways:01100 


15 


0.0001 


Allograft rejection:05330 


4 


0.0001 


MAPK signaling pathway:04010 


7 


0.0002 


Autoimmune thyroid disease:05320 


4 


0.0003 


VEGF signaling pathway:04370 


4 


0.0004 


Intestinal immune network for IgA 


3 


0.0010 


production:04672 






C21-Steroid hormone 


2 


0.0014 


metabolism:00140 






Graft-vs. host disease:05332 


3 


0.0022 


WIKI ENRICHMENT ANALYSIS 


Diurnally regulated genes with 


3 


0.0006 


circadian orthologs:WP1306 






B cell 


4 


0.0040 


receptor signaling pathway:WP285 






Cytokine and inflammatory 


2 


0.0040 


response:WP271 






Tryptophan metabolism:WP270 


2 


0.0128 


Regulation of actin 


3 


0.0201 


cytoskeleton:WP351 






TGF beta signaling pathway:WP505 


2 


0.0195 


MAPK signaling pathway:WP358 


3 


0.0238 


Kit receptor signaling pathway:WP147 


2 


0.0309 


IL-5 signaling pathway:WP44 


2 


0.0318 


GPCRs, class A rhodopsin-like:WP473 


3 


0.0564 



The effect of stress on differential expression of genes in the DRN 
is different depending on physical activity status 

The effect of stress on gene expression in the DRN is differ- 
ent depending on physical activity status. Figure 3 shows the 
number of genes changed by p < 0.05 that were differentially 
expressed between sedentary rats and physically active rats imme- 
diately following or 2 h post stress. Of these genes, 1028 were 
differentially expressed immediately following stress. Differential 
expression of 637 genes was observed 2 h after stress. Differential 
expression of 169 genes was observed in sedentary rats com- 
pared to physically active rats immediately following stress and 
differential expression of these genes was also present 2 h after 
stress. 

Genes that were differentially expressed at both time points 
(n = 169) were subjected to KEGG pathway analysis (Table 3) in 
order to identify functionally enriched pathway categories related 
to these genes. Stress differentially impacted the expression 
of genes related to functional categories including olfactory 
transduction, metabolic pathways, MAPK signaling pathways, 
and VEGF signaling pathways in physically active compared to 
sedentary rats. 
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Genes that were differentially expressed at both time points 
(n = 169) were also subjected to Wiki pathway analysis (Table 3). 
Stress differentially impacted the expression of genes related to 
functional categories including diurnally regulated genes with 
circadian orthologs, tryptophan metabolism, and GPCRs in 
physically active compared to sedentary rats. Various immune- 
associated pathway categories were also identified such as 
cytokine and inflammatory response pathways as well as signaling 
pathways for B-cell receptor, TGF-|i, and IL-5. 

The effect of wheel running and/or exposure to stress on the 
expression of select genes related to functionally enriched wiki 
pathway categories 

Select genes from Wiki functional pathway categories were tar- 
geted for analysis by AN OVA to reveal the main effect of wheel 
running, main effect of stress, and exercise x stress interaction on 
gene expression in the DRN immediately following and 2 h post 
stress. Table 4 shows the statistical results of the ANOVAs for each 
gene. Figure 4 displays the graphs for each AN OVA. Genes were 
selected from the following functional categories: diurnally reg- 
ulated genes with circadian orthologs, tryptophan metabolism, 
and various immune-related pathways including B-cell recep- 
tor signaling, TGF-|3 signaling, IL-5 signaling, and cytokines and 
inflammatory response. Within the diurnally regulated genes 
with circadian orthologs category, genes selected for ANOVA 
were Kruppel-like factor 9 (Klf9), protein phosphate 1 (Ppplr3c), 
and G0/G1 switch 2 (G0s2). Within the tryptophan metabolism 
category, genes selected for ANOVA were aldehyde dehydroge- 
nase 1 family (Aldhla2) and tryptophan 2,3-dioxygenase (Tdo2). 
Within the immune-related categories, genes selected for ANOVA 
included hematopoietic cell specific Lyn substrate 1 (Hclsl), 
phospholipase C, gamma 2 {Plcg2), Cas-Br-M ecotropic retro- 
viral transforming sequence b (Cblb), hemopoietic cell kinase 
(Hck), interleukin 4 (IL-4), transforming growth factor, beta 
1 (TGF-fil), and BMP and activin membrane-bound inhibitor 
(Bambi). 



WEIGHTED GENE CORRELATIONAL NETWORK ANALYSIS 

Hierarchical clustering of rats based on gene expression profiles 

In order to construct modules of highly correlated genes that 
were related to exercise and/or stress exposure, a WGCNA was 
performed. Hierarchical clustering was used to categorize rats 
based on their individual expression profile of the genes within 
the transcriptome (17,170). First, a dendogram was constructed 
to cluster rats based on gene expression profile. Rats that were 
closer in distance within the dendogram were considered to have 
a more closely related gene expression profile. After the dendo- 
gram was constructed, clusters of the dendogram were related 
to physical traits whereby experimental condition (exercise and 
stress) were considered physical traits. With the exception of the 
outliers (« = 7), rats fell into two main categories, home cage 
non-stressed rats on the left and rats exposed to stress on the 
right. All groups within the stress cluster contained both seden- 
tary and physically active rats with no visible grouping pattern 
by physical activity status. Within the home cage non-stressed 
control cluster, one stressed outlier was identified. Rats with the 
same physical activity status were clustered together within the 
non-stressed cluster. 

Identification of modules of coexpressed genes in the DRN 
correlated to 6 weeks of wheel running and/or exposure to stress 

Modules of highly coexpressed genes were identified that were 
also highly correlated (either positively or negatively) to exercise 
and/or stress. Eleven modules were derived from the transcrip- 
tome and were related to the various physical traits. Each module 
was assigned an arbitrary color. The number of genes con- 
tributing to each module was: yellow-199, blue-1077, purple-36, 
magenta-46, turquoise-3350, red-99, black-67, brown-373, green- 
153, pink-53, and grey-1 1,717. A correlation value and p-value 
associated with the correlational strength for each module-trait 
relationship was calculated. Modules of interest were identified 
based on statistical significance of the correlational strength (p < 
0.001 ). A more stringent cutoff for statistical significance was used 



Table 4 | Effect of exercise on stress-induced alterations in the expression of select genes in the DRN. 



FUNCTIONAL CATEGORY 



Gene Exercise Stress Exercise x stress 



DIURNALLY REGULATED W/CIRCADIAN ORTHOLOGS 



Km 


Fo 


38) 


= 0.012, p = 


0.915 


F(2, 38) 


= 7.842, p = 0.001 


F(2, 


38) 


= 4.552, p = 0.016 


Ppp1r3c 


Fa 


38) 


= 0.504, p = 


0.482 


F(2, 38) 


= 6.551, p = 0.003 


Fa. 


38) 


= 4.014, p = 0.026 


G0s2 


Fa 


38) 


= 3.239, p = 


0.079 


F(2, 38) 


= 3.813, p = 0.031 


F,2. 


38) 


= 4.408, p = 0.019 


TRYPTOPHAN METABOLISM 


Aldh1a2 


Fa 


38) 


= 2.713, p = 


0.107 


F(2, 38) 


= 0.384, p = 0.683 


Fa. 


38) 


= 2.891, p = 0.067 


Tdo2 


Fa 


38) 


= 0.658, p = 


0.422 


F(2, 38) 


= 1.243, p = 0.300 


F,2. 


38) 


= 5.022, p = 0.011 


IMMUNE-RELATED 


Hclsl 


Fa 


38) 


= 3.417, p = 


0.072 


F(2, 38) 


= 0.494, p = 0.614 


F,2. 


38) 


= 2.896, p = 0.067 


Plcg2 


Fa 


38) 


= 1 .898, p = 


0.176 


F(2, 38) 


= 0.129, p = 0.879 


F,2. 


38) 


= 4.352, p = 0.019 


Cblb 


Fa 


38) 


= 7.493, p = 


0.009 


F(2, 38) 


= 0.968, p = 0.389 


F,z 


38) 


= 3.179, p = 0.052 


Hck 


Fa 


38) 


= 8.404, p = 


0.006 


F(2, 38) 


= 1.858, p = 0.169 


Fa. 


38) 


= 4.010, p = 0.026 


IL-4 


Fa 


38) 


= 0.284, p = 


0.597 


F(2, 38) 


= 0.121, p = 0.886 


Fa. 


38) 


= 2.983, p = 0.062 


TGF-&1 


Fa 


38) 


= 0.528, p = 


0.472 


F(2, 38) 


= 22.91, p < 0.0001 


Fa. 


38) 


= 6.430, p = 0.003 


Bambi 


Fa 


38) 


= 8.404, p = 


0.006 


F(2, 38) 


= 1.858, p = 0.169 


Fa. 


38) 


= 4.010, p = 0.026 
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FIGURE 4 | Effect of exercise on stress-induced alterations in the expression of select genes in the DRN, organized by Wiki functional category. 
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given the wide range of p-values (p = 3e — 39 to p = 1.0). The 
grey module was excluded from analysis due to the large num- 
ber of genes contributing to the module, and therefore, lack of 
specificity. 

Of the 11 modules, 2 modules, the brown and black, were 
responsive to stress in both physically active and sedentary rats. 
The brown module was highly positively correlated to stress, indi- 
cating a strong increase in expression of genes within the brown 
module in response to stress. The physically active rats had a 
greater correlation value immediately following stress (0.98) com- 
pared to sedentary rats (0.83), suggesting that there was a more 
coordinated response among genes in the brown module in phys- 
ically active rats. The black module was also highly positively 
correlated to stress, indicating a strong increase in expression of 
genes within the black module in response to stress. For both time 
points, physically active rats had a greater correlation value (0.96 
and 0.83) compared to sedentary rats (0.8 and 0.72), suggesting 
that there was a more coordinated response among genes in the 
black module in physically active rats in response to stress. 

Five additional modules were also identified. Genes within 
the blue, purple, and green modules were responsive to stress 
only in the sedentary rats. The blue module was negatively corre- 
lated (—0.5) with stress in sedentary rats, indicating a decrease in 
expression of genes within the blue module in response to stress. 
The purple module was negatively correlated (—0.5 and —0.65) 
with stress in sedentary rats immediately following and 2 h post 
stress, indicating a decrease in expression of genes within the pur- 
ple module in response to stress. The green module was positively 
correlated (0.48) immediately following stress in sedentary rats, 
indicating an increase in expression of genes within the green 
module in response to stress. 

Expression of genes within the purple and turquoise mod- 
ules was associated with physical activity. The purple module 
was negatively correlated with physical activity (—0.48) indicat- 
ing a decrease in expression of genes within the purple module in 
response to wheel running. The turquoise module was positively 
correlated with physical activity (0.53) indicating an increase in 
expression of genes within the turquoise module in response to 
wheel running. 

Finally, the magenta module was positively correlated (0.53) 
with physically active rats 2 h post stress compared to sedentary 
rats 2 h post stress. This suggests that relative to sedentary rats, 
physically active rats have increased expression of genes within 
the magenta module 2 h following stress exposure. 

The effect of wheel running and/or exposure to stress on modules 
of coexpressed genes in the DRN 

Modules of interest (correlational strength p < 0.001) were also 
subjected to analysis by ANOVA in order to determine the effect 
of exercise and/or stress on alterations in the coexpression of 
genes in the DRN. Figure 5 shows the graphs of the ANOVA anal- 
ysis for each module. There were no statistically significant (p < 
0.05) effects of exercise, stress, or an exercise by stress interaction 
in the blue, purple, or green module. 

ANOVA analysis of the magenta module revealed a statistically 
significant main effect of exercise [-F(i, 38) = 6.588; p = 0.0143], 
but no significant main effect of stress or interaction. ANOVA 



analysis of the turquoise module also revealed a statistically sig- 
nificant main effect of exercise [i 7 (i_ 38) = 5.495; p = 0.0244], but 
no significant main effect of stress or interaction. ANOVA analysis 
of the black module module revealed a statistically significant 
main effect of stress [f(2, 38) = 56.872; p < 0.0001] and signifi- 
cant interaction [Fp, 28) = 3.431; p = 0.0427], but no significant 
main effect of exercise. ANOVA analysis of the brown module 
revealed a statistically significant main effect of stress [Fp, 38) = 
168.838; p < 0.001], but no significant main effect of exercise or 
interaction. 

KEGG functionally enriched pathway categories related to genes 
differentially coexpressed in the DRN following 6 weeks of wheel 
running and/or exposure to stress 

For modules found to be statistically significant through ANOVA 
analysis, genes that contributed to each module were imported 
into Web Gestaldt bioinformatics system and subjected to analysis 
with KEGG functional terms. Table 5 shows the top ten function- 
ally enriched pathway categories, for each module, related to the 
genes differentially coexpressed following 6 weeks of wheel run- 
ning and/or exposure to stress. Gene count refers to the number 
of genes from the data sets that contribute to each functional cat- 
egory. The p-value represents the statistical significance of each 
functionally enriched category identified. 

KEGG analysis revealed that genes within the black module 
were related to functional categories including prion diseases, 
Wnt signaling, chemokine signaling, toll-like receptor signaling, 
VEGF signaling, gonadotropin-releasing hormone (GnRH) sig- 
naling, apoptosis, and small cell lung cancer. Genes within the 
brown module were related to functional categories including 
Jak-Stat signaling, adipocytokine signaling, pathways in cancer, 
hematopoietic cell lineage, p53 signaling, focal adhesion, chronic 
myeloid leukemia, and ErbB signaling. Both the black and brown 
modules contained genes that were related to functional cate- 
gories of pathways involving cytokine-cytokine receptor interac- 
tion and MAPK signaling. 

Genes within the magenta and turquoise modules were related 
to functional categories including metabolic pathways and ubiq- 
uitin mediated proteolysis. Refer to Table 5 for other functional 
categories related to genes within these modules. 

Wiki functionally enriched pathway categories related to genes 
differentially coexpressed in the DRN following 6 weeks of wheel 
running and/or exposure to stress 

For modules found to be statistically significant through ANOVA 
analysis, genes that contributed to each module were imported 
into Web Gestaldt bioinformatics system and subjected to analysis 
with Wiki functional terms. Table 6 shows the top ten function- 
ally enriched pathway categories, for each module, related to the 
genes differentially coexpressed following 6 weeks of wheel run- 
ning and/or exposure to stress. Gene count refers to the number 
of genes from the data sets that contribute to each functional cat- 
egory. The p-value represents the statistical significance of each 
functionally enriched category identified. 

Wiki analysis revealed that genes within the brown module 
were related to functional categories including adipogenesis, IL- 
6 signaling, triacylglyceride synthesis, ErbB signaling, p38 MAPK 
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FIGURE 5 | Effect of exercise and/or stress on modules of coexpressed genes in the DRN. 
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Table 5 | KEGG functionally enriched pathway categories generated Table 6 | Wiki functionally enriched pathway categories generated 
from modules of genes correlated to exercise and/or stress. from modules of genes correlated to exercise and/or stress. 

Module Gene count P-value Module Gene count P-value 



FUNCTIONALLY ENRICHED PATHWAY 



MAGENTA MODULE 



Focal adhesion:04510 


3 


0.0014 


Type II diabetes mellitus:04930 


2 


0.0015 


Pancreatic cancer:05212 


2 


0.0031 


Colorectal cancer:05210 


2 


0.0047 


Pathways in cancer:05200 


3 


0.0061 


Neurotropin signaling pathway:04722 


2 


0.0099 


Insulin signaling pathway:04910 


2 


0.0099 


Ubiquitin mediated proteolysis:04120 


2 


0.0091 


Metabolic Pathways:01100 


2 


0.3704 




Metabolic Pathways:01100 


231 


3.92e-38 


Ribosome:03010 


43 


7.08e-24 


Ubiquitin mediated proteolysis:04120 


45 


9.21 e-18 


Axon guidance:04360 


42 


6.30e-15 


Huntington's disease:05016 


56 


2.63e-13 


Alzheimer's disease:05010 


57 


7.96e-13 


Oxidative phosphorylation:00190 


43 


3.01 e-12 


Regulation of actin cytoskeleton:04810 


51 


3.72e-12 


Spliceosome:03040 


37 


3.89e-12 


Cell cycle:04110 


37 


5.05e-12 




Cytokine-cytokine receptor interaction:04060 


4 


0.0005 


Prion Diseases:05020 


2 


0.0020 


Wnt signaling pathway:04310 


3 


0.0023 


MAPK signaling pathway:04010 


4 


0.0014 


Chemokine signaling pathway:04062 


3 


0.0039 


Toll-like receptor signaling pathway:04620 


2 


0.0127 


VEGF signaling pathway:04370 


2 


0.0087 


GnRH signaling pathway:04912 


2 


0.0132 


Apoptosis:04210 


2 


0.0137 


Small cell lung cancer:05222 


2 


0.0127 




MAPK signaling pathway:04010 


23 


3.94e-15 


Cytokine-cytokine receptor interaction:04060 


13 


9.53e-08 


Jak-Stat signaling pathway:04630 


10 


2.02e-06 


Adipocytokine signaling pathway:04920 


7 


4.19e-06 


Pathways in cancer:05200 


14 


6.26e-06 


Hematopoietic cell lineage:04640 


7 


1.07e-05 


P53 signaling pathway:04115 


6 


9.48e-05 


Focal adhesion:04510 


9 


0.0002 


Chronic myeloid leukemia:05220 


6 


0.0002 


ErbB signaling pathway:04012 


6 


0.0002 



signaling, and Wnt signaling and pluripotency. Genes within 
the black module were related to functional categories including 
hypertrophy model, small ligand GPCRs, prostaglandin synthesis 
and regulation, myometrial relaxation and contraction, diurnally 
regulated genes with circadian orthologs, and peptide GPCRs. 
Both the black and brown modules contained genes that were 



FUNCTIONALLY ENRICHED PATHWAY 





Insulin signaling:WP439 


3 


0.0005 


Apoptosis:WP1290 


2 


0.0039 


IL-3 


2 


0.0048 


signaling pathway:WP319 






Apoptosis mechanisms:WP284 


2 


0.0038 


TURQUOISE MODULE 


mRNA processing:WP529 


41 


2.45e-17 


Electron transport chain:WP59 


35 


6.06e-17 


TNF alpha NF-KB signaling pathway:WP457 


47 


2.77e-14 


EGFR1 signaling pathway:WP5 


42 


6.39e-12 


Regulation of actin cytoskeleton:WP351 


37 


2.25e-11 


TGF beta receptor signaling 


35 


1.04e-10 


pathway:WP362 






B cell 


38 


1.96e-10 


receptor signaling pathway:WP285 






G protein signaling pathway:WP73 


27 


2.60e-10 


Oxidative phosphorylation:WP1283 


19 


3.28e-09 


Proteasome degradation :WP302 


19 


1.51 e-08 


BLACK MODULE 


Hypertrophy Model:WP442 


3 


4.38e-06 


Insulin signaling:WP429 


4 


0.0001 


GPCRs, class A rhodopsin-like:WP473 


4 


0.0005 


Small ligand GPCRs:WP161 


2 


0.0004 


Prostaglandin synthesis and 


2 


0.0012 


regulation:WP303 






Myometrial relaxation and contraction 


3 


0.0018 


pathways:WP140 






MAPK signaling pathway:WP358 


3 


0.0022 


TGF beta receptor signaling 


3 


0.0016 


pathway:WP362 






Diurnally regulated genes with circadian 


2 


0.0022 


orthologs:WP1306 






Peptide GPCRs:WP131 


2 


0.0043 




MAPK signaling pathway:WP358 


15 


1.29e-11 


Adipogenesis:WP1 55 


12 


1.41 e-09 


Insulin signaling:WP429 


9 


1.04e-05 


TGF beta receptor signaling 


8 


3.95e-05 


pathway:WP362 






IL-6 signaling pathway:WP135 


7 


4.22e-05 


Triacylglyercide synthesis:WP356 


4 


4.69e-05 


ErbB signaling pathway:WP1299 


5 


7.64e-05 


P38 MAPK signaling pathway:WP294 


4 


0.0002 


GPCRs, class A rhodopsin-like:WP473 


9 


0.0002 


Wnt signaling pathway and 


6 


0.0002 


pluripotency:WP1288 








TNF alpha NF-KB signaling pathway:WP457 


5 


7.31 e-05 


Electron transport chain:WP59 


3 


0.0013 



/Continued) 



Frontiers in Behavioral Neuroscience 



www.frontiersin.org 



May 2013 | Volume 7 | Article 37 | 15 



Loughridge et al. 



Genetic networks of stress resistance 



Table 6 | Continued 



It/I 

IVIodulG 


Gene count 


r-value 


FUNCTIONALLY ENRICHED PATHWAY 




Androgen receptor signaling 


3 


0.0029 


pathway:WP68 






Cytoplasmic ribosomal proteins:WP30 


3 


0.0028 


Oxidative phosphorylation :WP1 283 


2 


0.0070 


Proteasome degradation :WP302 


2 


0.0081 


G1 to S cell cycle control :WP348 


2 


0.0110 


Wnt signaling pathway:WP375 


2 


0.0227 



related to functional categories of pathways in insulin signaling, 
MAPK signaling, GPCRs of class A rhodopsin-like, and TGF-f} 
receptor signaling. 

For the magenta and turquoise modules, functional categories 
associated with immune pathways such as IL-3 signaling, insulin 
signaling, TGF-f5 receptor signaling, B-cell receptor signaling, and 
TNF-a-NF-icB signaling were identified. 

DISCUSSION 
OVERALL THEMES 

The mechanism by which exercise protects against the behav- 
ioral consequences of inescapable stress is unknown. The current 
data suggest that rats with 6 weeks of prior access to a running 
wheel have a different physiological response to stress, as mea- 
sured by gene expression in the DRN, than sedentary rats. Here 
we report that (1) relative to home cage non-stressed controls, 
physically active rats have a greater number of genes differentially 
expressed in response to stress both immediately following and 
2 h after stress exposure than sedentary rats (2) modules made up 
of genes that are highly coexpressed and responsive to stress oper- 
ate in a more strongly coordinated manner in response to stress in 
physically active rats compared to sedentary rats (3) many of the 
stress-responsive genes within the DRN are known to be involved 
in various immune-related pathways, such as cytokine signaling 
and inflammatory processes. 

These data demonstrate that in response to stress, physically 
active rats mount a more active response, at the level of mRNA 
transcription in the DRN. Relative to home cage non-stressed 
controls, physically active rats had a greater number of genes 
altered by a LFC > ±1.1 and a greater number of genes sig- 
nificantly differentially altered by p < 0.05 than sedentary rats 
in response to stress. This is interesting because it suggests that 
the protective effect of exercise is not through a dampening of 
the stress response. Rather, physically active rats may mount 
a more robust response that functions, in concert, to protect 
the brain against the negative behavioral consequences of stress. 
Furthermore, physically active rats may respond to stress in a 
more efficient manner. This is demonstrated by the observation 
that modules made up of coexpressed genes that are highly stress- 
responsive, are more strongly upregulated in physically active rats 
in response to stress than sedentary rats. Given that these mod- 
ules are highly upregulated in response to stress despite physical 
activity status, it is possible that the genes within these modules 
are expressed in order to protect the DRN from stress-induced 



damage. Physically active rats have a more strongly coordi- 
nated coexpression of these stress-responsive, "protective" genes 
and this more effective synchronization may be important for 
exercise-induced stress resistance. 

Interestingly, many genes within the DRN that are altered in 
response to stress are involved in immune-related signaling pro- 
cesses including the signaling of proteins (MAPK) involved in the 
stimulation of proinflammatory factors, immune cell receptors 
(toll-like, B cell, T cell), cytokines (IL-3, IL-4, IL-5, IL-6, TGF-p\ 
TNF-ct) and cytokine receptors (TGF-fU), chemokines, and reg- 
ulatory pathways involved in the immune response to infection 
(NF-kB). These various immune-related functional categories 
were identified in both differential expression and WGCNA anal- 
ysis of genes and in both KEGG and Wild pathway databases. It is 
important to note that identification of specific functional path- 
ways does not necessarily imply that such processes are occurring 
within the DRN in response to exercise and/or stress. Rather, 
pathway identification is a means of organizing the thousands of 
genes that are significantly differentially expressed or coexpressed 
in response to exercise and/or stress, by functional relationships. 
Genes known to be involved in B cell receptor signaling, for 
example, were significantly altered in the DRN in response to 
stress. However, B cells are only present at very low levels in 
a healthy brain (Anthony et al., 2003). Thus, these genes are 
likely performing other functions. Additionally, it is likely that 
the immune-associated functional categories of genes within the 
KEGG and Wiki pathway databases were generated from data 
based on peripheral immune processes, and therefore, may not 
be an accurate representation of immune functions within the 
brain. However, given the evidence that physical activity has anti- 
inflammatory effects peripherally that confer protection against 
chronic disease (Gleeson et al., 201 1), elucidation of the nature of 
the stress-induced inflammatory response in the DRN of physi- 
cally active rats is of prime importance. It is possible that exercise 
confers protection by dampening the stress-induced inflamma- 
tory cascade in the brain. 

Overall, these data suggest that at the level of mRNA transcrip- 
tion in the DRN, there is a colossal response to stress. A history 
of physical activity, changes, but does not necessarily dampen this 
response. Furthermore, there is evidence of stress-induced induc- 
tion of inflammatory processes, though it is not clear whether 
there is an overall pro-inflammatory, anti-inflammatory, or bal- 
anced response. Additionally, it is not apparent if the overall 
inflammatory response is different depending on physical activ- 
ity status. Regardless, these data provide evidence for a stress- 
induced inflammatory response originating in a region of the 
brain implicated in stress-related mood disorders, and exem- 
plify the importance of investigating novel theories, such as the 
cytokine-induced hypothesis of depression. 

NOVEL TARGETS OF EXERCISE-INDUCED STRESS RESISTANCE 

The goal of this experiment was to identify novel gene targets 
of exercise-induced stress resistance. The "novel" component was 
considered to be of particular importance and therefore, the data 
were analyzed without the guidance of an a priori hypothesis. 
Differential gene expression analysis was employed to narrow 
the transcriptome of 17,170 genes to those most likely involved 
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in stress resistance. Specifically, contrasts were made between 
experimental groups to identify genes that were significantly 
differentially expressed by p < 0.05. The contrasts considered 
to be the most important for revealing exercise-induced stress 
resistance were the contrasts that revealed genes differentially 
expressed due to an exercise by stress interaction. More specifi- 
cally, we were interested in identifying genes that were differen- 
tially expressed at p < 0.05 in response to stress depending on 
the physical activity status of the rat. The contrasts that pro- 
vided this information were [(Sedstresso v. Sednc) v. (Runstresso v. 
Run H c)] and [(Sed St ress2 v. Sed H c) v. (Run St ress2 v. Run H c)]- The 
first contrast generated a list of 1028 genes that were differentially 
expressed immediately following stress in physically active rats 
compared to sedentary rats. The second contrast generated a list 
of 637 genes that were differentially expressed 2 h following stress 
in physically active rats compared to sedentary rats. Given the 
large number of genes that were altered, an additional constraint 
was placed on the data. Only those genes that were differentially 
expressed at both time points were considered (n = 169). We rea- 
soned that mRNA present immediately following stress and 2 h 
following stress was more likely to be translated into protein. The 
169 genes that were differentially expressed in response to stress 
in physically active rats compared to sedentary rats were subjected 
to both KEGG and Wiki analysis, in order to organize the list 
of genes by functional categories. It is important to note that of 
the 169 genes, only 69 genes were fitted to KEGG specific func- 
tional categories, and only 26 genes were fitted to Wiki specific 
functional categories. Thus, it is possible that important genes 
were lost due to the limitations inherent in classifying genes into 
categories. From an a priori perspective, visual inspection of the 
genes not assigned to a KEGG or Wiki functional category did not 
reveal any compelling genes, per se. We chose to focus on genes 
classified by the Wiki database because the functional categories 
seemed to be tailored to more specific processes (TGF-p" signaling 



Table 7 | Functional role of proteins encoded by genes of interest. 
FUNCTIONAL CATEGORY 



Gene Protein product: function 



DIURNALLY REGULATED W/CIRCADIAN ORTHOLOGS 

Klf9 Transcription factor: can inhibit or activate transcription 

Ppp1r3c Enzyme: participates in variety of cellular processes by reversible protein phosphorylation 

G0s2 Not clear: potential oncogene and regulator of latent HIV 

Aldh1a2 Enzyme: catalyzes the synthesis of retinoic acid from retinaldehyde 

Tdo2 Enzyme: catalyzes 1st and rate-limiting step in a major pathway of tryptophan metabolism, Lrtryptophan > 

n-formyl kynurenine 

Hclsl Substrate: role in antigen receptor signaling, potential role in regulation of gene expression 

Plcg2 Enzyme: catalyze hydrolysis of phospholipids 

Cblb Ligase: negatively regulates T-cell and B-cell receptors 

Hck Enzyme: may help couple Fc receptor to activation of respiratory burst, potential role in neutrophil migration 

and degranulation of neutrophils 

IL-4 Cytokine: forms a gene cluster w/ IL-3, IL-5, IM3 on chromosome 5q, costimulator of DNA synthesis, 

induces expression of MHC II on B-cells 
TGF-fi 1 Cytokine: controls proliferation, differentiation, regulation of other growth factors 

Bambi Receptor: related to type 1 receptors of TGF-p family 



and tryptophan metabolism pathways) compared to KEGG iden- 
tified processes (metabolic pathways and pathways in cancer). 
From the Wiki categories, 12 genes were selected that were related 
to functional categories including diurnally regulated genes, tryp- 
tophan metabolism, and inflammatory-related processes. These 
categories were specifically chosen because the circadian system, 
serotonergic circuits, and inflammation have all been implicated 
in having a role in mood disorders. It should be noted that of 
the 12 genes selected from Wiki functional categories, 6 of these 
genes were also present in the KEGG classification of functional 
categories. 

ANOVA analysis was performed in order to assess the dif- 
ferential effect of stress-induced changes within these genes in 
physically active compared to sedentary rats. Table 7 provides a 
brief description of the function associated with the protein that 
each gene encodes (Safran et al., 2010). 

Of particular interest are Tdo2 and TGF-fil (Figure 4). 
Tryptophan is the precursor to 5-HT and therefore, serves as 
a regulator of 5-HT synthesis. The Tdo2 gene encodes a major 
enzyme, tryptophan 2,3 dioxygenase (TDO), involved in trypto- 
phan metabolism. TDO degrades tryptophan along the kynure- 
nine pathway (for review, see Efimov et al., 2011). Interestingly, 
in rats exposed to stress, Tdo2 mRNA expression in the DRN is 
differentially altered depending on physical activity status [ExS, 
F(2, 38) = 5.022; p = 0.0116]. Fisher's PLSD post-hoc analyses 
revealed that exercise decreased baseline levels of Tdo2 and there 
was no effect of stress on Tdo2 expression in physically active 
rats. In contrast, sedentary rats had decreased Tdo2 expres- 
sion immediately following and 2h post stress. Stress -induced 
decreases in Tdo2 may lead to a reduction in TDO levels and 
decreased ability of rats to degrade tryptophan. This could con- 
tribute to the increased levels of 5-HT in the DRN following 
inescapable stress. Wheel running blocks stress-induced decreases 
in Tdo2 mRNA expression, and thus, may prevent stress-induced 
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decreases in TDO and the reduction of tryptophan metabolism in 
the DRN. 

TGF-fil was an additional gene of interest. The TGF-$1 gene 
encodes for the cytokine, TGF-pi. TGF-pi is a regulator of T cells, 
and is typically considered to have anti-inflammatory properties. 
TGF-P 1 has also been investigated in the context of depres- 
sion. Kim et al. (2007) found increased levels of TGF-fU in the 
plasma of patient's with major depression disorder and 6 weeks 
of treatment with antidepressants significantly decreased levels of 
TGF-pi. Our data suggests that TGF-fil mRNA in the DRN is par- 
ticularly sensitive to stress (p < 0.0001) and differentially altered 
depending on physical activity status [ExS, Fp, 38) = 6.430; p = 
0.0039]. Fisher's PLSD post-hoc analysis revealed that exercise 
increased baseline levels of TGF-fil mRNA. There was a decrease 
in TGF-fil mRNA expression in physically active rats 2h post 
stress. In contrast, sedentary rats had a large increase in TGF-$1 
immediately following stress. The significance of the impedance 
effect of exercise on stress-induced increases in TGF-$1 mRNA 
levels in the DRN is unclear. Sedentary rats may have a more pro- 
inflammatory response to stress, and therefore, the increase in 
TGF-$1 mRNA may occur in order to subdue pro-inflammatory 
responses. On the other hand, 5-HT treatment has been shown 
to increase TGF-fil mRNA expression in mesangial cells (Grewal 
et al., 1999) and therefore, increased TGF-$1 mRNA expression 
in DRN cells may be due to stress-induced increases in 5-HT lev- 
els in sedentary rats. However, it is important to consider that 
the functional role of TGF-P 1 maybe different in the brain com- 
pared to peripheral tissue. Given the association of TGF-P 1 with 
depression and antidepressant treatment, further examination of 
TGF-P 1 is warranted. 

It is important to consider that in the context of this particu- 
lar analysis, the identification of novel targets of exercise-induced 
stress resistance was restricted to those genes that were differ- 
entially expressed due to the interaction between exercise and 
stress. However, it is possible that the mechanism by which exer- 
cise confers protection is not opposite of the mechanism by 
which stress produces negative consequences. Exercise-induced 
stress resistance could occur through a non-stress responsive 
route. Furthermore, identification of stress-resistant genes relied 
solely on the genes being differentially expressed in physically 
active rats compared to sedentary rats exposed to stress. However, 
genes often work in coordinated manners to carry out physio- 
logical functions. In the absence of absolute differences in gene 
expression, differences in the coexpression of gene networks in 
response to stress may underlie the protective effect of exercise. 
Identification of hub genes critically important to the modules 
detected with WGCNA will address this possibility, and may lead 
to the identification of additional targets. (Identification of hub 
genes is discussed in greater detail in Future Directions.) 

ADVANTAGES AND LIMITATIONS OF MICR0ARRAY ANALYSIS 

The greatest advantage of microarray analysis is that it enables 
the simultaneous exploration of the expression of thousands of 
genes. Therefore, it is particularly useful in studying complex 
processes, such as the stress response, whereby thousands of 
genes are affected. When used in conjunction with laser cap- 
ture microarray technology, microarray has the potential to yield 



whole genome expression data about an organism's response to 
an environmental manipulation, such as exposure to stress or 
voluntary exercise, in a specific region of the brain or cell type 
within that region. It is important to point out that microar- 
ray analysis only provides information at the level of mRNA 
transcription, which is not necessarily indicative of protein pro- 
duction. Nevertheless, transcription initiation is the most widely 
used means of gene regulation in eukaryotes (Cox et al., 2012), 
and using an exploratory approach, microarray experiments serve 
as an important starting screen for the identification of novel tar- 
gets that can be analyzed in greater detail with other techniques 
(Coppola, 2011). 

The exploratory approach to microarray data analysis, how- 
ever, is not without its limitations. This approach relies solely 
on statistical significance to identify often thousands of genes 
that in turn, must be organized in such a way that allows 
for interpretation. The organization of statistically significant 
genes, usually by functional categories derived from bioinfor- 
matics databases, may fail to identify genes that play a crucial 
role in the regulation of a given process. The SHTi^receptor, 
for example, is thought to be critically involved in the mecha- 
nism by which inescapable stress produces sensitization of DRN 
5-HT neurons (Rozeske et al., 2011). Using in-situ hybridiza- 
tion, our lab has previously shown that wheel running increases 
5-HTu mRNA expression in the DRN (Greenwood et al, 2003, 
2005). An exploratory analysis of the current dataset did not iden- 
tify 5-HTiyi as a target of interest. However, using an a priori 
guided approach, ANOVA analysis of 5-HTia mRNA in the DRN 
revealed significant main effects of exercise [F(i i 38) = 6.250; p = 
0.0169] and stress [F (2 , 38 ) = 3.538; p = 0.0390], but no signif- 
icant interaction (Figure 6). Our data are consistent with the 
previous findings that 6 weeks of wheel running increases 5- 
HTla mRNA expression in the DRN (Greenwood et al., 2003). 
We also replicated our previous report that exposure to stress 
decreases 5-HTia mRNA expression in the DRN (Greenwood and 
Fleshner, 2011). 

Overall, both the advantages and limitations of microarray 
technology are the product of the colossal dataset that a microar- 
ray experiment yields. Organization of the data is required so 
that statistically significant differences observed in thousands of 
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genes can be focused into more manageable and interpretable 
gene sets. This process of data distillation, however, is not without 
a price and in the process, information may be lost. Given this 
reality, it is of paramount importance that multiple statistical 
and analytical strategies are executed. The traditional analysis of 
differential expression should be used in addition to the more 
sophisticated analysis of coexpression with WGCNA. Although 
both approaches return long lists of genes that must be fur- 
ther mined, each method assesses different regulatory processes. 
Finally, when identifying functional categories related to genes of 
interest, multiple bioinformatics databases should be explored. 
Inconsistencies in the information returned by bioinformatics 
databases (Soh et al., 2010) are pervasive, and using more than 
one database may provide more comprehensive results. 

FUTURE DIRECTIONS 

The overwhelming amount of data obtained in microarray exper- 
iments can be a challenge to manage, however, the results are 
powerful and provide researchers with a wealth of novel processes 
and genes to explore. To identify additional targets of exercise- 
induced stress resistance, hub genes should be identified within 
the modules of highly coexpressed genes detected by the WGCNA. 
Research suggests that a gene's position within a given network of 
genes, or module, is indicative of its functional significance to that 
module (Miller et al., 2008). More centralized genes, termed hub 
genes, are more likely to have critical roles in cellular function 
than peripheral genes (Miller et al., 2008). Therefore, centralized 
hub genes within the brown module may play a significant role in 
the effect of stress, while centralized hub genes within the black 
module (which was differentially altered in response to stress 
depending on physical activity status), may play a significant role 
in the protective effect of exercise against stress. Identification of 
hub genes within the brown and black modules could lead to the 
identification of additional therapeutic targets. 

Further analysis of significantly regulated genes with assign- 
ment of cell type may also add richly to the dataset. More specif- 
ically, a database containing information on the genes enriched 
in a given cell type, such as neurons, astrocytes, and oligoden- 
drocytes, can be used to assign cell type specificity to significantly 
differentially expressed or coexpressed genes (Cahoy et al, 2008). 
Given there was consistent evidence for inflammatory-related 
genes being upregulated in response to stress in the DRN, it may 
be insightful to identify the specific cell-type and source of this 
inflammatory-related mRNA expression. 

Future experiments should include validating the novel tar- 
gets identified with microarray analysis with PCR and/or in-situ 
hybridization. In-situ hybridization, in particular, could provide 
information on the anatomical specificity of the observed differ- 
ences in mRNA expression in response to exercise and/or stress 
within the DRN. Additionally, novel gene targets of exercise- 
induced stress resistance should be tested by means of behavioral 
pharmacology. That is manipulation of the proteins encoded by 
these genes with specific agonists or antagonists in the context of 
inescapable stress exposure and/or exercise may reveal informa- 
tion on the therapeutic potential of these genes. Delivery of phar- 
macological agonsits of targets of upregulated "stress-resistance 
genes" to sedentary rats, for example, would be expected to confer 



protection against stress-induced behaviors in these rats in the 
absence of exercise. 

Androgens and circadian regulationare additional topics that 
may warrant further investigation for having a role in the mecha- 
nisms by which exercise produces stress resistance. Our data sug- 
gest that pathways of genes related to androgen receptor signaling 
and diurnal regulation are differentially expressed in the DRN 
in physically active compared to sedentary rats following stress. 
Androgens promote neurogenesis (Spiritzer and Galea, 2007) and 
enhance cognitive function (Edinger and Frye, 2007). A recent 
study reported that mild exercise increased androgen synthesis in 
the hippocampus (Okamoto et al., 2012). It is possible that wheel 
running also increases androgen synthesis in the DRN and pro- 
vides stress-resistance through an androgen-mediated pathway. 
Additionally, disruption of the circadian rhythm has been impli- 
cated in mood disorders. Interestingly, stimulation of the DRN 
triggers the release of 5-HT by the suprachiasmatic nucleus (SCN) 
(Dudley et al., 1999), the brain region responsible for the gener- 
ation of circadian rhythms (Rusak and Zucker, 1979), and DRN 
stimulation induces circadian phase-resetting (Glass et al., 2000). 
Constraint of DRN 5-HT neurons in physically active rats could 
block stress-induced alterations in SCN output and provide stress 
resistance through the prevention of circadian rhythm disruption. 

CONCLUSIONS 

In conclusion, when the data are organized effectively, microarray 
experiments have the ability to yield a rich amount of information 
on the molecular activities underlying physiological processes. 
When used in combination with laser capture microdissection, 
this information can be obtained from a specific region or 
cell type within an organism. Thus microarray technology is 
particularly useful in studying the neurobiological mechanisms 
underlying the complex pathophysiology of stress-related mood 
disorders. This experiment was designed to reveal novel targets 
by which exercise produces resistance to stress-related mood dis- 
orders, specifically within the DRN, using microarray and laser 
capture microdissection technology. The current data reveal evi- 
dence for different profiles of gene expression in the DRN of 
physically active rats exposed to stress compared to sedentary 
rats exposed to stress. Physically active rats have a more active 
and more strongly coordinated response to stress than seden- 
tary rats. Specifically, Tdo2, a gene encoding an enzyme involved 
in tryptophan metabolism, may have a role in the mechanism 
by which exercise protects against the behavioral consequences 
of inescapable stress. In addition, an inflammatory-related gene 
encoding for the cytokine TGF-fil was particularly responsive 
to stress and this response was different depending on physical 
activity status. Overall, an inflammatory theme was revealed con- 
sistently across multiple analyses, suggesting a large effect of stress 
on inflammatory-related processes in cells of the DRN. The con- 
sequence of stress-induced inflammatory processes in the DRN 
should be further investigated. 
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